Loading libraries

Loading Data

Plots for Gross Photosynthesis - Cabral:

# mutate(ammonia_umolL = ifelse(ammonia_umolL = "<0.02", 0.01, ammonia_umolL)) %>% 

###################
## Silicate
####################

### confirm linear model 
GP_Cabral_silicate_linearmodel <- lmer(umol.cm2.hr ~ log(silicate_umolL+1) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Cabral") %>% 
  filter(P_R=="GP") %>% 
  filter(!is.infinite(umol.cm2.hr)))

anova(GP_Cabral_silicate_linearmodel)
## Type III Analysis of Variance Table with Satterthwaite's method
##                           Sum Sq  Mean Sq NumDF  DenDF F value  Pr(>F)  
## log(silicate_umolL + 1) 0.056267 0.056267     1 46.348  4.0643 0.04962 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
### confirm nonlinear model 
GP_Cabral_silicate_nonlinearmodel <- lmer(umol.cm2.hr ~ poly(log(silicate_umolL+1),2) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Cabral") %>% 
  filter(P_R=="GP") %>% 
  filter(!is.infinite(umol.cm2.hr)))

anova(GP_Cabral_silicate_nonlinearmodel)
## Type III Analysis of Variance Table with Satterthwaite's method
##                                   Sum Sq  Mean Sq NumDF  DenDF F value  Pr(>F)
## poly(log(silicate_umolL + 1), 2) 0.10844 0.054218     2 46.553  4.1984 0.02107
##                                   
## poly(log(silicate_umolL + 1), 2) *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#### both are significant so can add either trend line 

GP_Cabral_Silicate <- RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Cabral") %>% 
  filter(P_R=="GP") %>% 
  filter(!is.infinite(umol.cm2.hr)) %>% 
  ggplot(aes(x=silicate_umolL, 
             y=umol.cm2.hr)) + 
  geom_point(size=4, aes(color=as.factor(new_colonynumber))) +
  scale_color_manual(values=pnw_palette("Bay", n=9), 
                     guide="none") +
  theme_bw() + 
  coord_trans(x ="log") +
  theme(strip.background = element_rect(fill = "white")) +
 geom_smooth(method="lm", formula = "y~log(x)", color="black") +
 # geom_smooth(method="lm", formula = "y~poly(log(x),2)", color="black") +
 # geom_hline(yintercept = 0) +
  labs(x = "log silicate",
       y = "Rate (umol O2 cm2 hr-1)",
       title = "GP Rates for Cabral with Silicate") +
   theme(axis.text.x=element_text(size=18), 
        axis.text.y=element_text(size=18), 
        axis.title.x=element_text(size=18),
        axis.title.y=element_text(size=18), 
        plot.title=element_text(hjust=0.5, ## centers title 
                                size=20))
GP_Cabral_Silicate

###################
## SGD dils
####################

### confirm a linear model 
GP_Cabral_SGD_linearmodel <- lmer(umol.cm2.hr ~ log(SGD_number+1) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Cabral") %>% 
  filter(P_R=="GP") %>% 
  filter(!is.infinite(umol.cm2.hr)))

anova(GP_Cabral_SGD_linearmodel)
## Type III Analysis of Variance Table with Satterthwaite's method
##                       Sum Sq  Mean Sq NumDF  DenDF F value  Pr(>F)  
## log(SGD_number + 1) 0.045966 0.045966     1 46.338   3.266 0.07722 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
### confirm a nonlinear model 
GP_Cabral_SGD_nonlinearmodel <- lmer(umol.cm2.hr ~ poly(log(SGD_number+1),2) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Cabral") %>% 
  filter(P_R=="GP") %>% 
  filter(!is.infinite(umol.cm2.hr)))

anova(GP_Cabral_SGD_nonlinearmodel)
## Type III Analysis of Variance Table with Satterthwaite's method
##                               Sum Sq  Mean Sq NumDF  DenDF F value Pr(>F)  
## poly(log(SGD_number + 1), 2) 0.10382 0.051908     2 45.272  3.9695 0.0258 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
GP_Cabral_SGDDils <- RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Cabral") %>% 
  filter(P_R=="GP") %>% 
  filter(!is.infinite(umol.cm2.hr)) %>% 
  ggplot(aes(x=SGD_number+1, 
             y=umol.cm2.hr)) + 
  geom_point(size=4, aes(color=as.factor(new_colonynumber))) +
  scale_color_manual(values=pnw_palette("Bay", n=9), 
                     guide="none") +
  theme_bw() + 
  coord_trans(x ="log") +
  theme(strip.background = element_rect(fill = "white")) +
  geom_smooth(method="lm", formula = "y~poly(log(x),2)", color="black") +
 # geom_hline(yintercept = 0) +
  labs(x = "log SGD Dils",
       y = "Rate (umol O2 cm2 hr-1)",
       title = "GP Rates for Cabral with SGD Dils") +
   theme(axis.text.x=element_text(size=18), 
        axis.text.y=element_text(size=18), 
        axis.title.x=element_text(size=18),
        axis.title.y=element_text(size=18), 
        plot.title=element_text(hjust=0.5, ## centers title 
                                size=20))
GP_Cabral_SGDDils

###################
## SGD dils WITHOUT POTENTIAL OUTLIERS
####################
GP_Cabral_SGDDils_nooutlier <- RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Cabral") %>% 
  filter(sample_ID!= "Cabral_Col6_Dil9_Light") %>% 
  filter(P_R=="GP") %>% 
  filter(!is.infinite(umol.cm2.hr)) %>% 
  ggplot(aes(x=SGD_number+1, 
             y=umol.cm2.hr)) + 
  geom_point(size=4, aes(color=as.factor(new_colonynumber))) +
 scale_color_manual(values=pnw_palette("Bay", n=9)) +
  theme_bw() + 
  coord_trans(x ="log") +
  theme(strip.background = element_rect(fill = "white")) +
  #geom_smooth(method="lm", formula = "y~poly(log(x),2)", color="black") +
  geom_smooth(method="lm", formula = "y~log(x)", color="black") +
 # geom_hline(yintercept = 0) +
  labs(x = "log SGD Dils",
       y = "Rate (umol O2 cm2 hr-1)",
       title = "GP Rates for Cabral with SGD Dils") +
   theme(axis.text.x=element_text(size=18), 
        axis.text.y=element_text(size=18), 
        axis.title.x=element_text(size=18),
        axis.title.y=element_text(size=18), 
        plot.title=element_text(hjust=0.5, ## centers title 
                                size=20))
GP_Cabral_SGDDils_nooutlier

### confirm a linear model 
GP_Cabral_SGD_linearmodel_nooutlier <- lmer(umol.cm2.hr ~ log(SGD_number+1) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Cabral") %>% 
  filter(sample_ID!= "Cabral_Col6_Dil9_Light") %>% 
  filter(P_R=="GP") %>% 
  filter(!is.infinite(umol.cm2.hr)))

anova(GP_Cabral_SGD_linearmodel_nooutlier)
## Type III Analysis of Variance Table with Satterthwaite's method
##                       Sum Sq  Mean Sq NumDF DenDF F value  Pr(>F)  
## log(SGD_number + 1) 0.081847 0.081847     1    50  6.2821 0.01549 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
### confirm a nonlinear model 
GP_Cabral_SGD_nonlinearmodel_nooutlier <- lmer(umol.cm2.hr ~ poly(log(SGD_number+1),2) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Cabral") %>% 
  filter(P_R=="GP") %>% 
  filter(sample_ID!= "Cabral_Col6_Dil9_Light") %>% 
  filter(!is.infinite(umol.cm2.hr)))

anova(GP_Cabral_SGD_nonlinearmodel_nooutlier)
## Type III Analysis of Variance Table with Satterthwaite's method
##                               Sum Sq  Mean Sq NumDF DenDF F value  Pr(>F)  
## poly(log(SGD_number + 1), 2) 0.11464 0.057321     2    49  4.5402 0.01553 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Plots for Gross Photosynthesis - Varari

###################
## Silicate
####################

### confirm linear model 
GP_Varari_silicate_linearmodel <- lmer(umol.cm2.hr ~ log(silicate_umolL+1) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Varari") %>% 
  filter(P_R=="GP") %>%
  filter(!is.infinite(umol.cm2.hr)))

anova(GP_Varari_silicate_linearmodel)
## Type III Analysis of Variance Table with Satterthwaite's method
##                           Sum Sq  Mean Sq NumDF DenDF F value Pr(>F)
## log(silicate_umolL + 1) 0.015463 0.015463     1 62.55  0.7003 0.4059
### confirm nonlinear model 
GP_Varari_silicate_nonlinearmodel <- lmer(umol.cm2.hr ~ poly(log(silicate_umolL+1),2) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Varari") %>% 
  filter(P_R=="GP") %>% 
  filter(!is.infinite(umol.cm2.hr)))

anova(GP_Varari_silicate_nonlinearmodel)
## Type III Analysis of Variance Table with Satterthwaite's method
##                                   Sum Sq  Mean Sq NumDF  DenDF F value  Pr(>F)
## poly(log(silicate_umolL + 1), 2) 0.15527 0.077633     2 61.577   3.925 0.02487
##                                   
## poly(log(silicate_umolL + 1), 2) *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#### ONLY significant in nonlinear model 

GP_Varari_Silicate <- RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Varari") %>% 
  filter(P_R=="GP") %>% 
  filter(!is.infinite(umol.cm2.hr)) %>% 
  ggplot(aes(x=silicate_umolL, 
             y=umol.cm2.hr)) + 
  geom_point(size=4, aes(color=as.factor(new_colonynumber))) +
  scale_color_manual(values=pnw_palette("Bay", n=9), 
                     guide="none") +
  theme_bw() + 
  coord_trans(x ="log") +
  theme(strip.background = element_rect(fill = "white")) +
 # geom_smooth(method="lm", formula = "y~log(x)", color="black") +
  geom_smooth(method="lm", formula = "y~poly(log(x),2)", color="black") +
 # geom_hline(yintercept = 0) +
  labs(x = "log silicate",
       y = "Rate (umol O2 cm2 hr-1)",
       title = "GP Rates for Varari with Silicate") +
   theme(axis.text.x=element_text(size=18), 
        axis.text.y=element_text(size=18), 
        axis.title.x=element_text(size=18),
        axis.title.y=element_text(size=18), 
        plot.title=element_text(hjust=0.5, ## centers title 
                                size=20))
GP_Varari_Silicate

###################
## SGD dils
####################

### confirm a linear model 
GP_Varari_SGD_linearmodel <- lmer(umol.cm2.hr ~ log(SGD_number+1) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Varari") %>% 
  filter(P_R=="GP") %>% 
  filter(!is.infinite(umol.cm2.hr)))

anova(GP_Varari_SGD_linearmodel)
## Type III Analysis of Variance Table with Satterthwaite's method
##                       Sum Sq  Mean Sq NumDF  DenDF F value Pr(>F)
## log(SGD_number + 1) 0.010297 0.010297     1 60.065  0.4631 0.4988
### confirm a nonlinear model 
GP_Varari_SGD_nonlinearmodel <- lmer(umol.cm2.hr ~ poly(log(SGD_number+1),2) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Varari") %>% 
  filter(P_R=="GP") %>% 
  filter(!is.infinite(umol.cm2.hr)))

anova(GP_Varari_SGD_nonlinearmodel)
## Type III Analysis of Variance Table with Satterthwaite's method
##                                Sum Sq  Mean Sq NumDF  DenDF F value Pr(>F)
## poly(log(SGD_number + 1), 2) 0.089681 0.044841     2 59.102  2.1119   0.13
GP_Varari_SGDDils <- RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Varari") %>% 
  filter(P_R=="GP") %>% 
  filter(!is.infinite(umol.cm2.hr)) %>% 
  ggplot(aes(x=SGD_number+1, 
             y=umol.cm2.hr)) + 
  geom_point(size=4, aes(color=as.factor(new_colonynumber))) +
  scale_color_manual(values=pnw_palette("Bay", n=9), 
                     guide="none") +
  theme_bw() + 
  coord_trans(x ="log") +
  theme(strip.background = element_rect(fill = "white")) +
#  geom_smooth(method="lm", formula = "y~poly(log(x),2)", color="black") +
 # geom_hline(yintercept = 0) +
  labs(x = "log SGD Dils",
       y = "Rate (umol O2 cm2 hr-1)",
       title = "GP Rates for Varari with SGD Dils") +
   theme(axis.text.x=element_text(size=18), 
        axis.text.y=element_text(size=18), 
        axis.title.x=element_text(size=18),
        axis.title.y=element_text(size=18), 
        plot.title=element_text(hjust=0.5, ## centers title 
                                size=20))
GP_Varari_SGDDils

Plots for Respiration - Cabral:

###################
## SGD dils
####################

### confirm a linear model 
R_Cabral_SGD_linearmodel <- lmer(umol.cm2.hr ~ log(SGD_number+1) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Cabral") %>% 
  filter(P_R=="R") %>% 
  filter(!is.infinite(umol.cm2.hr)))

anova(R_Cabral_SGD_linearmodel)
## Type III Analysis of Variance Table with Satterthwaite's method
##                        Sum Sq   Mean Sq NumDF  DenDF F value Pr(>F)
## log(SGD_number + 1) 9.037e-08 9.037e-08     1 46.085       0 0.9962
### confirm a nonlinear model 
R_Cabral_SGD_nonlinearmodel <- lmer(umol.cm2.hr ~ poly(log(SGD_number+1),2) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Cabral") %>% 
  filter(P_R=="R") %>% 
  filter(!is.infinite(umol.cm2.hr)))

anova(R_Cabral_SGD_nonlinearmodel)
## Type III Analysis of Variance Table with Satterthwaite's method
##                                Sum Sq   Mean Sq NumDF  DenDF F value Pr(>F)
## poly(log(SGD_number + 1), 2) 0.011554 0.0057771     2 45.076  1.5607 0.2211
R_Cabral_SGDDils <- RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Cabral") %>% 
  filter(P_R=="R") %>% 
  filter(!is.infinite(umol.cm2.hr)) %>% 
  ggplot(aes(x=SGD_number+1, 
             y=umol.cm2.hr)) + 
  geom_point(size=4, aes(color=as.factor(new_colonynumber))) +
  scale_color_manual(values=pnw_palette("Bay", n=9), 
                     guide="none") +
  theme_bw() + 
  coord_trans(x ="log") +
  theme(strip.background = element_rect(fill = "white")) +
 # geom_smooth(method="lm", formula = "y~poly(log(x),2)", color="black") +
 # geom_hline(yintercept = 0) +
  labs(x = "log SGD Dils",
       y = "Rate (umol O2 cm2 hr-1)",
       title = "R Rates for Cabral with SGD Dils") +
   theme(axis.text.x=element_text(size=18), 
        axis.text.y=element_text(size=18), 
        axis.title.x=element_text(size=18),
        axis.title.y=element_text(size=18), 
        plot.title=element_text(hjust=0.5, ## centers title 
                                size=20))
R_Cabral_SGDDils

Plots for Respiration - Varari

###################
## SGD dils
####################

### confirm a linear model 
R_Varari_SGD_linearmodel <- lmer(umol.cm2.hr ~ log(SGD_number+1) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Varari") %>% 
  filter(P_R=="R") %>% 
  filter(!is.infinite(umol.cm2.hr)))

anova(R_Varari_SGD_linearmodel)
## Type III Analysis of Variance Table with Satterthwaite's method
##                        Sum Sq   Mean Sq NumDF  DenDF F value Pr(>F)
## log(SGD_number + 1) 0.0066251 0.0066251     1 59.962  2.1962 0.1436
### confirm a nonlinear model 
R_Varari_SGD_nonlinearmodel <- lmer(umol.cm2.hr ~ poly(log(SGD_number+1),2) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Varari") %>% 
  filter(P_R=="R") %>% 
  filter(!is.infinite(umol.cm2.hr)))

anova(R_Varari_SGD_nonlinearmodel)
## Type III Analysis of Variance Table with Satterthwaite's method
##                                Sum Sq   Mean Sq NumDF  DenDF F value Pr(>F)
## poly(log(SGD_number + 1), 2) 0.012828 0.0064138     2 58.984  2.1691 0.1233
R_Varari_SGDDils <- RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Varari") %>% 
  filter(P_R=="R") %>% 
  filter(!is.infinite(umol.cm2.hr)) %>% 
  ggplot(aes(x=SGD_number+1, 
             y=umol.cm2.hr)) + 
  geom_point(size=4, aes(color=as.factor(new_colonynumber))) +
  scale_color_manual(values=pnw_palette("Bay", n=9), 
                     guide="none") +
  theme_bw() + 
  coord_trans(x ="log") +
  theme(strip.background = element_rect(fill = "white")) +
#  geom_smooth(method="lm", formula = "y~poly(log(x),2)", color="black") +
 # geom_hline(yintercept = 0) +
  labs(x = "log SGD Dils",
       y = "Rate (umol O2 cm2 hr-1)",
       title = "R Rates for Varari with SGD Dils") +
   theme(axis.text.x=element_text(size=18), 
        axis.text.y=element_text(size=18), 
        axis.title.x=element_text(size=18),
        axis.title.y=element_text(size=18), 
        plot.title=element_text(hjust=0.5, ## centers title 
                                size=20))
R_Varari_SGDDils

Plots for Calcification - Cabral:

###################
## pH
####################

### confirm nonlinear model 
C_Cabral_pH_nonlinearmodel <- lmer(umol.cm2.hr ~ poly(log(new_pH+1),2) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Cabral") %>% 
  filter(P_R=="C") %>% 
  filter(!is.infinite(umol.cm2.hr)))

anova(C_Cabral_pH_nonlinearmodel)
## Type III Analysis of Variance Table with Satterthwaite's method
##                            Sum Sq  Mean Sq NumDF  DenDF F value  Pr(>F)  
## poly(log(new_pH + 1), 2) 0.043887 0.021943     2 47.288  3.2971 0.04566 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
C_Cabral_pH <- RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Cabral") %>% 
  filter(P_R=="C") %>% 
  filter(!is.infinite(umol.cm2.hr)) %>% 
  ggplot(aes(x=new_pH, 
             y=umol.cm2.hr)) + 
  geom_point(size=4, aes(color=as.factor(new_colonynumber))) +
  scale_color_manual(values=pnw_palette("Bay", n=9), 
                     guide="none") +
  theme_bw() + 
  coord_trans(x ="log") +
  theme(strip.background = element_rect(fill = "white")) +
# geom_smooth(method="lm", formula = "y~log(x)", color="black") +
 geom_smooth(method="lm", formula = "y~poly(log(x),2)", color="black") +
 # geom_hline(yintercept = 0) +
  labs(x = "log pH",
       y = "Rate (umol O2 cm2 hr-1)",
       title = "C Rates for Cabral with pH") +
   theme(axis.text.x=element_text(size=18), 
        axis.text.y=element_text(size=18), 
        axis.title.x=element_text(size=18),
        axis.title.y=element_text(size=18), 
        plot.title=element_text(hjust=0.5, ## centers title 
                                size=20))
C_Cabral_pH

###################
## SGD dils
####################

### confirm a linear model 
C_Cabral_SGD_linearmodel <- lmer(umol.cm2.hr ~ log(SGD_number+1) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Cabral") %>% 
  filter(P_R=="C") %>% 
  filter(!is.infinite(umol.cm2.hr)))

anova(C_Cabral_SGD_linearmodel)
## Type III Analysis of Variance Table with Satterthwaite's method
##                       Sum Sq  Mean Sq NumDF  DenDF F value Pr(>F)
## log(SGD_number + 1) 0.015433 0.015433     1 46.137  2.1383 0.1504
### confirm a nonlinear model 
C_Cabral_SGD_nonlinearmodel <- lmer(umol.cm2.hr ~ poly(log(SGD_number+1),2) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Cabral") %>% 
  filter(P_R=="C") %>% 
  filter(!is.infinite(umol.cm2.hr)))

anova(C_Cabral_SGD_nonlinearmodel)
## Type III Analysis of Variance Table with Satterthwaite's method
##                               Sum Sq   Mean Sq NumDF  DenDF F value Pr(>F)
## poly(log(SGD_number + 1), 2) 0.01618 0.0080902     2 45.126  1.0993 0.3418
C_Cabral_SGDDils <- RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Cabral") %>% 
  filter(P_R=="C") %>% 
  filter(!is.infinite(umol.cm2.hr)) %>% 
  ggplot(aes(x=SGD_number+1, 
             y=umol.cm2.hr)) + 
  geom_point(size=4, aes(color=as.factor(new_colonynumber))) +
  scale_color_manual(values=pnw_palette("Bay", n=9), 
                     guide="none") +
  theme_bw() + 
  coord_trans(x ="log") +
  theme(strip.background = element_rect(fill = "white")) +
  #geom_smooth(method="lm", formula = "y~poly(log(x),2)", color="black") +
 # geom_hline(yintercept = 0) +
  labs(x = "log SGD Dils",
       y = "Rate (umol O2 cm2 hr-1)",
       title = "C Rates for Cabral with SGD Dils") +
   theme(axis.text.x=element_text(size=18), 
        axis.text.y=element_text(size=18), 
        axis.title.x=element_text(size=18),
        axis.title.y=element_text(size=18), 
        plot.title=element_text(hjust=0.5, ## centers title 
                                size=20))
C_Cabral_SGDDils

Plots for Calcification - Varari

###################
## Phosphate
####################

### confirm nonlinear model 
C_Varari_phosphate_nonlinearmodel <- lmer(umol.cm2.hr ~ poly(log(phosphate_umolL+1),2) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Varari") %>% 
  filter(P_R=="C") %>% 
    filter(phosphate_umolL<1) %>% 
  filter(!is.infinite(umol.cm2.hr)))

anova(C_Varari_phosphate_nonlinearmodel)
## Type III Analysis of Variance Table with Satterthwaite's method
##                                    Sum Sq  Mean Sq NumDF  DenDF F value
## poly(log(phosphate_umolL + 1), 2) 0.16319 0.081594     2 60.327  7.2454
##                                     Pr(>F)   
## poly(log(phosphate_umolL + 1), 2) 0.001513 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
C_Varari_Phosphate <- RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Varari") %>% 
  filter(P_R=="C") %>% 
   filter(phosphate_umolL<1) %>% 
  filter(!is.infinite(umol.cm2.hr)) %>% 
  ggplot(aes(x=phosphate_umolL, 
             y=umol.cm2.hr)) + 
  geom_point(size=4, aes(color=as.factor(new_colonynumber))) +
  scale_color_manual(values=pnw_palette("Bay", n=9), 
                     guide="none") +
  theme_bw() + 
  coord_trans(x ="log") +
  theme(strip.background = element_rect(fill = "white")) +
 # geom_smooth(method="lm", formula = "y~log(x)", color="black") +
  geom_smooth(method="lm", formula = "y~poly(log(x),2)", color="black") +
 # geom_hline(yintercept = 0) +
  labs(x = "log phosphate",
       y = "Rate (umol O2 cm2 hr-1)",
       title = "C Rates for Varari with Phosphate") +
   theme(axis.text.x=element_text(size=18), 
        axis.text.y=element_text(size=18), 
        axis.title.x=element_text(size=18),
        axis.title.y=element_text(size=18), 
        plot.title=element_text(hjust=0.5, ## centers title 
                                size=20))
C_Varari_Phosphate

###################
## TA 
####################

### confirm a nonlinear model 
C_Varari_SGD_nonlinearmodel <- lmer(umol.cm2.hr ~ poly(log(TA_initial+1),2) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Varari") %>% 
    filter(TA_initial<2550) %>% #### IF CHOOSE TO REMOVE HIGH POINTS, ONLY SIGNIFICANT IN NONLINEAR (otherwise significant in linear too)
  filter(P_R=="C") %>% 
  filter(!is.infinite(umol.cm2.hr)))

anova(C_Varari_SGD_nonlinearmodel)
## Type III Analysis of Variance Table with Satterthwaite's method
##                               Sum Sq  Mean Sq NumDF  DenDF F value  Pr(>F)  
## poly(log(TA_initial + 1), 2) 0.10522 0.052612     2 59.782  3.6216 0.03276 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
C_Varari_TA <- RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Varari") %>% 
  filter(P_R=="C") %>% 
  # filter(TA_initial<2550) %>% 
  filter(!is.infinite(umol.cm2.hr)) %>% 
  ggplot(aes(x=TA_initial+1, 
             y=umol.cm2.hr)) + 
  geom_point(size=4, aes(color=as.factor(new_colonynumber))) +
  scale_color_manual(values=pnw_palette("Bay", n=9), 
                     guide="none") +
  theme_bw() + 
  coord_trans(x ="log") +
  theme(strip.background = element_rect(fill = "white")) +
 geom_smooth(method="lm", formula = "y~poly(log(x),2)", color="black") +
  geom_hline(yintercept = 0) +
  labs(x = "log TA",
       y = "Rate (umol O2 cm2 hr-1)",
       title = "C Rates for Varari with TA (Minus the potential outliers)") +
   theme(axis.text.x=element_text(size=18), 
        axis.text.y=element_text(size=18), 
        axis.title.x=element_text(size=18),
        axis.title.y=element_text(size=18), 
        plot.title=element_text(hjust=0.5, ## centers title 
                                size=20))
C_Varari_TA

###################
## SGD dils
####################

### confirm a linear model 
C_Varari_SGD_linearmodel <- lmer(umol.cm2.hr ~ log(SGD_number+1) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>%
  filter(site=="Varari") %>% 
  filter(P_R=="C") %>% 
  filter(!is.infinite(umol.cm2.hr)))

anova(C_Varari_SGD_linearmodel)
## Type III Analysis of Variance Table with Satterthwaite's method
##                      Sum Sq Mean Sq NumDF  DenDF F value   Pr(>F)   
## log(SGD_number + 1) 0.11101 0.11101     1 60.079  7.3097 0.008907 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
### confirm a nonlinear model 
C_Varari_SGD_nonlinearmodel <- lmer(umol.cm2.hr ~ poly(log(SGD_number+1),2) + (1|new_colonynumber), data=RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Varari") %>%
  filter(P_R=="C") %>% 
  filter(!is.infinite(umol.cm2.hr)))

anova(C_Varari_SGD_nonlinearmodel)
## Type III Analysis of Variance Table with Satterthwaite's method
##                               Sum Sq  Mean Sq NumDF  DenDF F value  Pr(>F)  
## poly(log(SGD_number + 1), 2) 0.12033 0.060163     2 59.118  3.9338 0.02489 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
C_Varari_SGDDils <- RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(site=="Varari") %>% 
  filter(P_R=="C") %>% 
  filter(!is.infinite(umol.cm2.hr)) %>% 
  ggplot(aes(x=SGD_number+1, 
             y=umol.cm2.hr)) + 
  geom_point(size=4, aes(color=as.factor(new_colonynumber))) +
  scale_color_manual(values=pnw_palette("Bay", n=9), 
                     guide="none") +
  theme_bw() + 
  coord_trans(x ="log") +
  theme(strip.background = element_rect(fill = "white")) +
geom_smooth(method="lm", formula = "y~log(x)", color="black") +
  geom_hline(yintercept = 0) +
  labs(x = "log SGD Dils",
       y = "Rate (umol O2 cm2 hr-1)",
       title = "NEC Rates for Varari with SGD Dils") +
   theme(axis.text.x=element_text(size=18), 
        axis.text.y=element_text(size=18), 
        axis.title.x=element_text(size=18),
        axis.title.y=element_text(size=18), 
        plot.title=element_text(hjust=0.5, ## centers title 
                                size=20))
C_Varari_SGDDils

Plot nutrient levels with SGD dils to see how different the nutrients actually are across the dilutions

Plot how many corals are dead vs alive

Net_CalcvsDissolution <- RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(P_R=="C") %>% 
  mutate(SGD_number = as.character(SGD_number), 
         rate_category = ifelse(umol.cm2.hr > 0, "positive", "negative")) %>% 
  filter(!is.infinite(umol.cm2.hr)) %>% 
  ggplot(aes(x=SGD_number, 
             y=umol.cm2.hr,
             fill=rate_category)) + 
  geom_col(width = 0.8,  colour = NA) + 
  #scale_fill_manual(values =  ifelse (umol.cm2.hr > 0, "red", "black")) +
  scale_fill_manual(values=c(positive='black',negative='red')) +
  theme_bw() + 
  geom_hline(yintercept = 0) +
  facet_wrap(~site) +
 #coord_trans(x ="log") +
  labs(x = "SGD values",
       y = "Rate (umol O2 cm2 hr-1)",
       title = "Net Calcifying vs Net Dissolving Corals per SGD Dilution") +
   theme(axis.text.x=element_text(size=14), 
        axis.text.y=element_text(size=18), 
        axis.title.x=element_text(size=18),
        axis.title.y=element_text(size=18), 
        plot.title=element_text(hjust=0.5, ## centers title 
                                size=20))
Net_CalcvsDissolution

Plot GP:R as a function of SGD

## First make new df
GP_RRatio_bySGDDils <- RespoR_Normalized_Full_withNEC_nutrients %>% 
  dplyr::select(!"sample_ID") %>% 
  filter(!is.infinite(umol.cm2.hr)) %>% 
  group_by(P_R, SGD_number, date, site, new_colonynumber) %>% 
  filter(P_R!="C", P_R!="NP") %>% 
  #summarize(GP_R = GP/R)
  pivot_wider(names_from = P_R, values_from = umol.cm2.hr) %>% 
  mutate(GP_R = GP/R) 

### plots 
GP_R_bySGD_bothsites <- GP_RRatio_bySGDDils %>% 
  #filter(site=="Varari") %>% 
  mutate(new_colonynumber= as.factor(new_colonynumber)) %>% 
  ggplot(aes(x=(SGD_number+1), 
             y=GP_R)) + 
  geom_point(size=3, aes(color=new_colonynumber)) +
  theme_bw() + 
  coord_trans(x="log") +
  facet_wrap(~site, scales = "free_y") +
  geom_smooth(method="lm", formula = "y~log(x)") +
  labs(x = "SGD Dilutions",
       y = "Gross Photosynthesis/Respiration Rates",
       title = "SGD by GP:R Relationship for Varari") +
   theme(axis.text.x=element_text(size=12), 
        axis.text.y=element_text(size=18), 
        axis.title.x=element_text(size=18),
        axis.title.y=element_text(size=18), 
        plot.title=element_text(hjust=0.5, ## centers title 
                                size=20))
GP_R_bySGD_bothsites

## model testing 
model_GP_R_SGD_Varari <- lmer(GP_R ~ log(SGD_number+1) + (1|new_colonynumber), data = GP_RRatio_bySGDDils %>%
                                        filter(site=="Varari"))
anova(model_GP_R_SGD_Varari)
## Type III Analysis of Variance Table with Satterthwaite's method
##                      Sum Sq Mean Sq NumDF  DenDF F value  Pr(>F)  
## log(SGD_number + 1) 0.26869 0.26869     1 60.003  3.6316 0.06148 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
model_GP_R_SGD_Cabral <- lmer(GP_R ~ log(SGD_number+1) + (1|new_colonynumber), data = GP_RRatio_bySGDDils %>%
                                        filter(site=="Cabral"))
anova(model_GP_R_SGD_Cabral)
## Type III Analysis of Variance Table with Satterthwaite's method
##                      Sum Sq Mean Sq NumDF  DenDF F value Pr(>F)  
## log(SGD_number + 1) 0.28731 0.28731     1 46.031  5.3534 0.0252 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Make correlation matrix for nutrients

## data 
corr_matrix_envparams <- RespoR_Normalized_Full_withNEC_nutrients %>% 
  filter(P_R=="R") %>% 
  mutate(Salinity=salinity) %>% 
  select(!c("salinity", "umol.cm2.hr", "SGD_number", "new_colonynumber", "sample_ID", "site", "date", "P_R", "ammonia_umolL")) ## removed ammonia for rn because have less than values and causing NAs

##visualize it 
ggcorr(corr_matrix_envparams, method = c("everything", "pearson")) 

###############################
### Dani B's methods #######
cor_mat <- round(cor(corr_matrix_envparams),4)
head(cor_mat)
##                 phosphate_umolL silicate_umolL NN_umolL  temp_C  new_pH
## phosphate_umolL          1.0000         0.2483   0.2789  0.1626 -0.0260
## silicate_umolL           0.2483         1.0000   0.1296 -0.1146 -0.1927
## NN_umolL                 0.2789         0.1296   1.0000  0.1277 -0.0638
## temp_C                   0.1626        -0.1146   0.1277  1.0000  0.2110
## new_pH                  -0.0260        -0.1927  -0.0638  0.2110  1.0000
## TA_initial               0.2644         0.4922   0.6274 -0.1274 -0.2502
##                 TA_initial Salinity
## phosphate_umolL     0.2644  -0.2635
## silicate_umolL      0.4922  -0.6914
## NN_umolL            0.6274   0.2955
## temp_C             -0.1274  -0.1727
## new_pH             -0.2502  -0.0230
## TA_initial          1.0000   0.0597
#melt the correlation matrix means it reassembles data frame to be more effective to complete corr matrix
#to long format
melted_cormat <- melt(cor_mat)
head(melted_cormat)
##              Var1            Var2   value
## 1 phosphate_umolL phosphate_umolL  1.0000
## 2  silicate_umolL phosphate_umolL  0.2483
## 3        NN_umolL phosphate_umolL  0.2789
## 4          temp_C phosphate_umolL  0.1626
## 5          new_pH phosphate_umolL -0.0260
## 6      TA_initial phosphate_umolL  0.2644
#visulaize the correlation matrix in general
ggplot(data = melted_cormat, aes(x=Var1, y=Var2, fill=value)) + 
  geom_tile()

# Get lower and upper triangle of the correlation matrix
#Note that, a correlation matrix has redundant information. We’ll use the functions below to set half of it to NA
get_lower_tri<-function(cor_mat){
  cormat[upper.tri(cor_mat)] <- NA
  return(cor_mat)
}
# Get upper triangle of the correlation matrix
get_upper_tri <- function(cor_mat){
  cor_mat[lower.tri(cor_mat)]<- NA
  return(cor_mat)
}

#apply upper tri calculation to graphc
upper_tri <- get_upper_tri(cor_mat)
upper_tri
##                 phosphate_umolL silicate_umolL NN_umolL  temp_C  new_pH
## phosphate_umolL               1         0.2483   0.2789  0.1626 -0.0260
## silicate_umolL               NA         1.0000   0.1296 -0.1146 -0.1927
## NN_umolL                     NA             NA   1.0000  0.1277 -0.0638
## temp_C                       NA             NA       NA  1.0000  0.2110
## new_pH                       NA             NA       NA      NA  1.0000
## TA_initial                   NA             NA       NA      NA      NA
## Salinity                     NA             NA       NA      NA      NA
##                 TA_initial Salinity
## phosphate_umolL     0.2644  -0.2635
## silicate_umolL      0.4922  -0.6914
## NN_umolL            0.6274   0.2955
## temp_C             -0.1274  -0.1727
## new_pH             -0.2502  -0.0230
## TA_initial          1.0000   0.0597
## Salinity                NA   1.0000
#melt the correlation matrix
#melt the correlation data and drop the rows with NA values 
melted_cormat <- melt(upper_tri, na.rm = TRUE)

#heatmap of correlation matrix
#negative correlations are in purple color and positive correlations in red
#scale_fill_gradient2 is used with the argument limit = c(-1,1) as correlation coefficients range from -1 to 1
#coord_fixed() : this function ensures that one unit on the x-axis is the same length as one unit on the y-axis
ggplot(data = melted_cormat, aes(Var2, Var1, fill = value))+
  geom_tile(color = "white")+
  scale_fill_gradient2(low = "midnightblue", high = "firebrick4", mid = "white", 
                       midpoint = 0, limit = c(-1,1), space = "Lab", 
                       name="Pearson\nCorrelation") +
  theme_minimal()+ 
  theme(axis.text.x = element_text(angle = 45, vjust = 1, 
                                   size = 12, hjust = 1))+
  coord_fixed()

melted_cormat <- melted_cormat %>% mutate_at(vars(starts_with("value")), funs(round(., 2)))

# Create a ggheatmap with basic characteristics, etc. 
ggheatmap <- ggplot(melted_cormat, aes(Var2, Var1, fill = value))+
  geom_tile(color = "white")+
  scale_fill_gradient2(low = "firebrick3", mid = "white", high = "dodgerblue3", 
                       midpoint = 0, limit = c(-1,1), space = "Lab", 
                       name="Pearson\nCorrelation") +
  theme_minimal()+ # minimal theme
  theme(axis.text.x = element_text(angle = 45, vjust = 1, 
                                   size = 12, hjust = 1))+
  coord_fixed()


# Print the heatmap
print(ggheatmap)

#add correlation coefficients to the heatmap
#geom_text() to add the correlation coefficients on the graph
#guides() to change the position of the legend title
#if else statement in melted data frame to quotes of black and white to adjust text color 
corr_matrix_SGD_params <- ggheatmap + 
  geom_text(aes(Var2, Var1, label = value), color = "black", size = 6) +
  theme(axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        axis.text.x = element_text(size = 18, face="bold", color="black"),
        axis.text.y = element_text(size = 18, face="bold", color="black"),
        panel.grid.major = element_blank(),
        panel.border = element_blank(),
        panel.background = element_blank(),
        axis.ticks = element_blank(),
        legend.justification = c(0.8, 0),
        legend.title = element_text(size = 18, face="bold", color="black"),
        legend.text = element_text(size = 20, face="bold", color="black"),
        legend.position = c(0.48, 0.75),
        legend.direction = "horizontal") +
  guides(fill = guide_colorbar(barwidth = 12, barheight = 2, 
                               title.position = "top", title.hjust = 0.5, title.vjust = 1.0))

corr_matrix_SGD_params

Results of back calculating